function Xrad = sepzeri(fun, xmin, xmax, nx, verb)

% Imposta numumero suddivisioni iniziali se non definite
if nargin < 5, verb = true; end
if nargin < 4, nx = 20; end

% Inizializzazione variabili d'uscita
Xrad = [];

% Inizializzazione parametri di ricerca
x = linspace(xmin, xmax, nx); % Vettori dei potenziali limiti delle radici
f = feval(fun, x); % Vettore delle f(x) nei punti degli intervalli

% Ricerca degli intervalli e loro visualizzazione
for k = 1 : length(x) - 1
    % Verifica se il segno di f(x) cambia
    if sign(f(k)) ~= sign(f(k + 1))
        % Salvataggio degli estremi dell'intervallo
        Xrad(end + 1, :) = x(k:k+1);
    end
end

% Stampa su warning se non trovato
if isempty(Xrad)
    warining('Nessuna separazione degli zeri trovata.');
end

% Visualizzazione grafica della soluzione trovata se richiesta
if verb
    % Disegna funzione con almeno 101 punti o 4 volte quelli iniziali
    xp = linspace(xmin, xmax, max(nx * 4 + 1, 101));
    yp = feval(fun, xp);
    
    % Box di visualizzazione grafica scalata (5% dell'originale)
    % per la visualizzazione di un rettangolo per ogni radice
    
    ytop = 0.05 * max(yp); ybot = 0.05 * min(yp); % Coordinate y
    ybox = [ybot ytop ytop ybot ybot]; % Box
    
    c = [0.5 0.5 0.5]; % Colore
    indx = [1 1 2 2 1]; % Indice dei punti Xrad fissato k della box
    
    % Disegno prima i rettangoli delle radici e poi la funzione
    figure;
    clf;
    hold on
    for i = 1 : size(Xrad, 1)
        fill(Xrad(i, indx), ybox, c);
    end
    hold off
    
    % Funzione
    hold on; h = plot(xp, yp, 'r-'); set(h, 'LineWidth', 2); hold off;
    grid on; xlabel('x'); ylabel('f(x)'); title('Separazione radici');
    axis tight;
end